home *** CD-ROM | disk | FTP | other *** search
/ STraTOS 1997 April & May / STraTOS 1 - 1997 April & May.iso / CD01 / INTERNET / SITES / LITTLE / P3SRC.ZIP / ATARI / HCMPLX.C < prev    next >
Encoding:
C/C++ Source or Header  |  1996-05-01  |  31.0 KB  |  1,523 lines
  1. /*****************************************************************************
  2. *                hcmplx.c
  3. *
  4. *  This module implements hypercomplex Julia fractals.
  5. *
  6. *  This file was written by Pascal Massimino.
  7. *
  8. *  from Persistence of Vision(tm) Ray Tracer
  9. *  Copyright 1996 Persistence of Vision Team
  10. *---------------------------------------------------------------------------
  11. *  NOTICE: This source code file is provided so that users may experiment
  12. *  with enhancements to POV-Ray and to port the software to platforms other
  13. *  than those supported by the POV-Ray Team.  There are strict rules under
  14. *  which you are permitted to use this file.  The rules are in the file
  15. *  named POVLEGAL.DOC which should be distributed with this file. If
  16. *  POVLEGAL.DOC is not available or for more info please contact the POV-Ray
  17. *  Team Coordinator by leaving a message in CompuServe's Graphics Developer's
  18. *  Forum.  The latest version of POV-Ray may be found there as well.
  19. *
  20. * This program is based on the popular DKB raytracer version 2.12.
  21. * DKBTrace was originally written by David K. Buck.
  22. * DKBTrace Ver 2.0-2.12 were written by David K. Buck & Aaron A. Collins.
  23. *
  24. *****************************************************************************/
  25.  
  26. #include "frame.h"
  27. #include "povray.h"
  28. #include "vector.h"
  29. #include "povproto.h"
  30. #include "fractal.h"
  31. #include "spheres.h"
  32. #include "hcmplx.h"
  33.  
  34.  
  35.  
  36. /*****************************************************************************
  37. * Local preprocessor defines
  38. ******************************************************************************/
  39.  
  40. #ifndef Fractal_Tolerance
  41. #define Fractal_Tolerance 1e-8
  42. #endif
  43.  
  44.  
  45. #define HMult(xr,yr,zr,wr,x1,y1,z1,w1,x2,y2,z2,w2)        \
  46.     (xr) = (x1) * (x2) - (y1) * (y2) - (z1) * (z2) + (w1) * (w2);   \
  47.     (yr) = (y1) * (x2) + (x1) * (y2) - (w1) * (z2) - (z1) * (w2);   \
  48.     (zr) = (z1) * (x2) - (w1) * (y2) + (x1) * (z2) - (y1) * (w2);   \
  49.     (wr) = (w1) * (x2) + (z1) * (y2) + (y1) * (z2) + (x1) * (w2);
  50.  
  51. #define HSqr(xr,yr,zr,wr,x,y,z,w)         \
  52.     (xr) = (x) * (x) - (y) * (y) - (z) * (z) + (w) * (w) ;  \
  53.     (yr) = 2.0 * ( (x) * (y) - (z) * (w) );     \
  54.     (zr) = 2.0 * ( (z) * (x) - (w) * (y) );       \
  55.     (wr) = 2.0 * ( (w) * (x) + (z) * (y) );
  56.  
  57.  
  58.  
  59. /*****************************************************************************
  60. * Local typedefs
  61. ******************************************************************************/
  62.  
  63.  
  64.  
  65. /*****************************************************************************
  66. * Static functions
  67. ******************************************************************************/
  68.  
  69. static int HReciprocal PARAMS((DBL * xr, DBL * yr, DBL * zr, DBL * wr, DBL x, DBL y, DBL z, DBL w));
  70. static void HFunc PARAMS((DBL * xr, DBL * yr, DBL * zr, DBL * wr, DBL x, DBL y, DBL z, DBL w, FRACTAL *f));
  71.  
  72.  
  73.  
  74. /*****************************************************************************
  75. * Local variables
  76. ******************************************************************************/
  77.  
  78. static CMPLX exponent = {0,0};
  79.  
  80. /******** Computations with Hypercomplexes **********/
  81.  
  82.  
  83.  
  84. /*****************************************************************************
  85. *
  86. * FUNCTION
  87. *
  88. * INPUT
  89. *
  90. * OUTPUT
  91. *
  92. * RETURNS
  93. *
  94. * AUTHOR
  95. *
  96. *   Pascal Massimino
  97. *
  98. * DESCRIPTION
  99. *
  100. * CHANGES
  101. *
  102. ******************************************************************************/
  103.  
  104. static int HReciprocal(xr, yr, zr, wr, x, y, z, w)
  105. DBL *xr, *yr, *zr, *wr, x, y, z, w;
  106. {
  107.   DBL det, mod, xt_minus_yz;
  108.  
  109.   det = ((x - w) * (x - w) + (y + z) * (y + z)) * ((x + w) * (x + w) + (y - z) * (y - z));
  110.  
  111.   if (det == 0.0)
  112.   {
  113.     return (-1);
  114.   }
  115.  
  116.   mod = (x * x + y * y + z * z + w * w);
  117.  
  118.   xt_minus_yz = x * w - y * z;
  119.  
  120.   *xr = (x * mod - 2 * w * xt_minus_yz) / det;
  121.   *yr = (-y * mod - 2 * z * xt_minus_yz) / det;
  122.   *zr = (-z * mod - 2 * y * xt_minus_yz) / det;
  123.   *wr = (w * mod - 2 * x * xt_minus_yz) / det;
  124.  
  125.   return (0);
  126. }
  127.  
  128.  
  129.  
  130. /*****************************************************************************
  131. *
  132. * FUNCTION Hfunc
  133. *
  134. * INPUT 4D Hypercomplex number, pointer to fractal object
  135. *
  136. * OUTPUT  calculates the 4D generalization of fractal->Complex_Function
  137. *
  138. * RETURNS void
  139. *
  140. * AUTHOR
  141. *
  142. *   Pascal Massimino
  143. *
  144. * DESCRIPTION
  145. *   Hypercomplex numbers allow generalization of any complex->complex
  146. *   function in a uniform way. This function implements a general
  147. *   unary 4D function based on the corresponding complex function. 
  148. *
  149. * CHANGES
  150. *  Generalized to use Complex_Function()   TW 
  151. *
  152. ******************************************************************************/
  153.  
  154. static void HFunc(xr, yr, zr, wr, x, y, z, w, f)
  155. DBL *xr, *yr, *zr, *wr, x, y, z, w;
  156. FRACTAL *f;
  157. {
  158.   CMPLX a, b, ra, rb;
  159.   
  160.   /* convert to duplex form */
  161.   a.x = x - w;
  162.   a.y = y + z;
  163.   b.x = x + w;
  164.   b.y = y - z;
  165.   
  166.   if(f->Sub_Type == PWR_STYPE)
  167.   {
  168.      exponent = f->exponent;
  169.   }
  170.   
  171.   /* apply function to each part */
  172.   Complex_Function(&ra, &a, f);
  173.   Complex_Function(&rb, &b, f);
  174.  
  175.   /* convert back */
  176.   *xr = .5 * (ra.x + rb.x);
  177.   *yr = .5 * (ra.y + rb.y);
  178.   *zr = .5 * (ra.y - rb.y);
  179.   *wr = .5 * (rb.x - ra.x);
  180. }
  181.  
  182.  
  183.  
  184. /*****************************************************************************
  185. *
  186. * FUNCTION
  187. *
  188. * INPUT
  189. *
  190. * OUTPUT
  191. *
  192. * RETURNS
  193. *
  194. * AUTHOR
  195. *
  196. *   Pascal Massimino
  197. *
  198. * DESCRIPTION
  199. *
  200. * CHANGES
  201. *
  202. ******************************************************************************/
  203.  
  204. /*------------------ z2 Iteration method ------------------*/
  205.  
  206. int Iteration_HCompl(IPoint, HCompl)
  207. VECTOR IPoint;
  208. FRACTAL *HCompl;
  209. {
  210.   int i;
  211.   DBL yz, xw;
  212.   DBL Exit_Value;
  213.   DBL x, y, z, w;
  214.  
  215.   x = Sx[0] = IPoint[X];
  216.   y = Sy[0] = IPoint[Y];
  217.   z = Sz[0] = IPoint[Z];
  218.   w = Sw[0] = (HCompl->SliceDist
  219.                   - HCompl->Slice[X]*x 
  220.                   - HCompl->Slice[Y]*y 
  221.                   - HCompl->Slice[Z]*z)/HCompl->Slice[T]; 
  222.   
  223.   Exit_Value = HCompl->Exit_Value;
  224.  
  225.   for (i = 1; i <= HCompl->n; ++i)
  226.   {
  227.     yz = y * y + z * z;
  228.     xw = x * x + w * w;
  229.  
  230.     if ((xw + yz) > Exit_Value)
  231.     {
  232.       return (FALSE);
  233.     }
  234.  
  235.     Sx[i] = xw - yz + HCompl->Julia_Parm[X];
  236.     Sy[i] = 2.0 * (x * y - z * w) + HCompl->Julia_Parm[Y];
  237.     Sz[i] = 2.0 * (x * z - w * y) + HCompl->Julia_Parm[Z];
  238.     Sw[i] = 2.0 * (x * w + y * z) + HCompl->Julia_Parm[T];
  239.  
  240.     w = Sw[i];
  241.     x = Sx[i];
  242.  
  243.     z = Sz[i];
  244.     y = Sy[i];
  245.   }
  246.  
  247.   return (TRUE);
  248. }
  249.  
  250.  
  251.  
  252. /*****************************************************************************
  253. *
  254. * FUNCTION
  255. *
  256. * INPUT
  257. *
  258. * OUTPUT
  259. *
  260. * RETURNS
  261. *
  262. * AUTHOR
  263. *
  264. *   Pascal Massimino
  265. *
  266. * DESCRIPTION
  267. *
  268. * CHANGES
  269. *
  270. ******************************************************************************/
  271.  
  272. int D_Iteration_HCompl(IPoint, HCompl, Dist)
  273. VECTOR IPoint;
  274. FRACTAL *HCompl;
  275. DBL *Dist;
  276. {
  277.   int i;
  278.   DBL yz, xw;
  279.   DBL Exit_Value, F_Value, Step;
  280.   DBL x, y, z, w;
  281.   VECTOR H_Normal;
  282.  
  283.   x = Sx[0] = IPoint[X];
  284.   y = Sy[0] = IPoint[Y];
  285.   z = Sz[0] = IPoint[Z];
  286.   w = Sw[0] = (HCompl->SliceDist
  287.                   - HCompl->Slice[X]*x 
  288.                   - HCompl->Slice[Y]*y 
  289.                   - HCompl->Slice[Z]*z)/HCompl->Slice[T]; 
  290.  
  291.   Exit_Value = HCompl->Exit_Value;
  292.  
  293.   for (i = 1; i <= HCompl->n; ++i)
  294.   {
  295.     yz = y * y + z * z;
  296.     xw = x * x + w * w;
  297.  
  298.     if ((F_Value = xw + yz) > Exit_Value)
  299.     {
  300.       Normal_Calc_HCompl(H_Normal, i - 1, HCompl);
  301.  
  302.       VDot(Step, H_Normal, Direction);
  303.  
  304.       if (Step < -Fractal_Tolerance)
  305.       {
  306.         Step = -2.0 * Step;
  307.  
  308.         if ((F_Value > Precision * Step) && (F_Value < 30 * Precision * Step))
  309.         {
  310.           *Dist = F_Value / Step;
  311.  
  312.           return (FALSE);
  313.         }
  314.       }
  315.  
  316.       *Dist = Precision;
  317.  
  318.       return (FALSE);
  319.     }
  320.  
  321.     Sx[i] = xw - yz + HCompl->Julia_Parm[X];
  322.     Sy[i] = 2.0 * (x * y - z * w) + HCompl->Julia_Parm[Y];
  323.     Sz[i] = 2.0 * (x * z - w * y) + HCompl->Julia_Parm[Z];
  324.     Sw[i] = 2.0 * (x * w + y * z) + HCompl->Julia_Parm[T];
  325.  
  326.     w = Sw[i];
  327.     x = Sx[i];
  328.  
  329.     z = Sz[i];
  330.     y = Sy[i];
  331.   }
  332.  
  333.   *Dist = Precision;
  334.  
  335.   return (TRUE);
  336. }
  337.  
  338.  
  339.  
  340. /*****************************************************************************
  341. *
  342. * FUNCTION
  343. *
  344. * INPUT
  345. *
  346. * OUTPUT
  347. *
  348. * RETURNS
  349. *
  350. * AUTHOR
  351. *
  352. *   Pascal Massimino
  353. *
  354. * DESCRIPTION
  355. *
  356. * CHANGES
  357. *
  358. ******************************************************************************/
  359.  
  360. void Normal_Calc_HCompl(Result, N_Max, fractal)
  361. VECTOR Result;
  362. int N_Max;
  363. FRACTAL *fractal;
  364. {
  365.   DBL n1, n2, n3, n4;
  366.   int i;
  367.   DBL x, y, z, w;
  368.   DBL xx, yy, zz, ww;
  369.   DBL Pow;
  370.  
  371.   /*
  372.    * Algebraic properties of hypercomplexes allows simplifications in
  373.    * computations...
  374.    */
  375.  
  376.   x = Sx[0];
  377.   y = Sy[0];
  378.   z = Sz[0];
  379.   w = Sw[0];
  380.  
  381.   Pow = 2.0;
  382.  
  383.   for (i = 1; i < N_Max; ++i)
  384.   {
  385.  
  386.     /*
  387.      * For a map z->f(z), f depending on c, one must perform here :
  388.      *
  389.      * (x,y,z,w) * df/dz(Sx[i],Sy[i],Sz[i],Sw[i]) -> (x,y,z,w) ,
  390.      *
  391.      * up to a constant.
  392.      */
  393.  
  394.     /******************* Case z->z^2+c *****************/
  395.  
  396.     /* the df/dz part needs no work */
  397.  
  398.     HMult(xx, yy, zz, ww, Sx[i], Sy[i], Sz[i], Sw[i], x, y, z, w);
  399.  
  400.     w = ww;
  401.     z = zz;
  402.     y = yy;
  403.     x = xx;
  404.  
  405.     Pow *= 2.0;
  406.   }
  407.  
  408.   n1 = Sx[N_Max] * Pow;
  409.   n2 = Sy[N_Max] * Pow;
  410.   n3 = Sz[N_Max] * Pow;
  411.   n4 = Sw[N_Max] * Pow;
  412.  
  413.   Result[X] = x * n1 + y * n2 + z * n3 + w * n4;
  414.   Result[Y] = -y * n1 + x * n2 - w * n3 + z * n4;
  415.   Result[Z] = -z * n1 - w * n2 + x * n3 + y * n4;
  416. }
  417.  
  418.  
  419.  
  420. /*****************************************************************************
  421. *
  422. * FUNCTION
  423. *
  424. * INPUT
  425. *
  426. * OUTPUT
  427. *
  428. * RETURNS
  429. *
  430. * AUTHOR
  431. *
  432. *   Pascal Massimino
  433. *
  434. * DESCRIPTION
  435. *
  436. * CHANGES
  437. *
  438. ******************************************************************************/
  439.  
  440. int F_Bound_HCompl(Ray, Fractal, Depth_Min, Depth_Max)
  441. RAY *Ray;
  442. FRACTAL *Fractal;
  443. DBL *Depth_Min, *Depth_Max;
  444. {
  445.   return (Intersect_Sphere(Ray, Fractal->Center, Fractal->Radius_Squared, Depth_Min, Depth_Max));
  446. }
  447.  
  448. /****************************************************************/
  449. /*--------------------------- z3 -------------------------------*/
  450. /****************************************************************/
  451.  
  452.  
  453.  
  454. /*****************************************************************************
  455. *
  456. * FUNCTION
  457. *
  458. * INPUT
  459. *
  460. * OUTPUT
  461. *
  462. * RETURNS
  463. *
  464. * AUTHOR
  465. *
  466. *   Pascal Massimino
  467. *
  468. * DESCRIPTION
  469. *
  470. * CHANGES
  471. *
  472. ******************************************************************************/
  473.  
  474. int Iteration_HCompl_z3(IPoint, HCompl)
  475. VECTOR IPoint;
  476. FRACTAL *HCompl;
  477. {
  478.   int i;
  479.   DBL Norm, xx, yy, zz, ww;
  480.   DBL Exit_Value;
  481.   DBL x, y, z, w;
  482.  
  483.   x = Sx[0] = IPoint[X];
  484.   y = Sy[0] = IPoint[Y];
  485.   z = Sz[0] = IPoint[Z];
  486.   w = Sw[0] = (HCompl->SliceDist
  487.                   - HCompl->Slice[X]*x 
  488.                   - HCompl->Slice[Y]*y 
  489.                   - HCompl->Slice[Z]*z)/HCompl->Slice[T]; 
  490.  
  491.   Exit_Value = HCompl->Exit_Value;
  492.  
  493.   for (i = 1; i <= HCompl->n; ++i)
  494.   {
  495.     Norm = x * x + y * y + z * z + w * w;
  496.  
  497.     /* is this test correct ? */
  498.     if (Norm > Exit_Value)
  499.     {
  500.       return (FALSE);
  501.     }
  502.  
  503.     /*************** Case: z->z^2+c *********************/
  504.     HSqr(xx, yy, zz, ww, x, y, z, w);
  505.  
  506.     x = Sx[i] = xx + HCompl->Julia_Parm[X];
  507.     y = Sy[i] = yy + HCompl->Julia_Parm[Y];
  508.     z = Sz[i] = zz + HCompl->Julia_Parm[Z];
  509.     w = Sw[i] = ww + HCompl->Julia_Parm[T];
  510.  
  511.   }
  512.  
  513.   return (TRUE);
  514. }
  515.  
  516.  
  517.  
  518. /*****************************************************************************
  519. *
  520. * FUNCTION
  521. *
  522. * INPUT
  523. *
  524. * OUTPUT
  525. *
  526. * RETURNS
  527. *
  528. * AUTHOR
  529. *
  530. *   Pascal Massimino
  531. *
  532. * DESCRIPTION
  533. *
  534. * CHANGES
  535. *
  536. ******************************************************************************/
  537.  
  538. int D_Iteration_HCompl_z3(IPoint, HCompl, Dist)
  539. VECTOR IPoint;
  540. FRACTAL *HCompl;
  541. DBL *Dist;
  542. {
  543.   int i;
  544.   DBL xx, yy, zz, ww;
  545.   DBL Exit_Value, F_Value, Step;
  546.   DBL x, y, z, w;
  547.   VECTOR H_Normal;
  548.  
  549.   x = Sx[0] = IPoint[X];
  550.   y = Sy[0] = IPoint[Y];
  551.   z = Sz[0] = IPoint[Z];
  552.   w = Sw[0] = (HCompl->SliceDist
  553.                   - HCompl->Slice[X]*x 
  554.                   - HCompl->Slice[Y]*y 
  555.                   - HCompl->Slice[Z]*z)/HCompl->Slice[T]; 
  556.  
  557.   Exit_Value = HCompl->Exit_Value;
  558.  
  559.   for (i = 1; i <= HCompl->n; ++i)
  560.   {
  561.     F_Value = x * x + y * y + z * z + w * w;
  562.  
  563.     if (F_Value > Exit_Value)
  564.     {
  565.       Normal_Calc_HCompl_z3(H_Normal, i - 1, HCompl);
  566.  
  567.       VDot(Step, H_Normal, Direction);
  568.  
  569.       if (Step < -Fractal_Tolerance)
  570.       {
  571.         Step = -2.0 * Step;
  572.  
  573.         if ((F_Value > Precision * Step) && (F_Value < 30 * Precision * Step))
  574.         {
  575.           *Dist = F_Value / Step;
  576.  
  577.           return (FALSE);
  578.         }
  579.       }
  580.  
  581.       *Dist = Precision;
  582.  
  583.       return (FALSE);
  584.     }
  585.  
  586.     /*************** Case: z->z^2+c *********************/
  587.  
  588.     HSqr(xx, yy, zz, ww, x, y, z, w);
  589.  
  590.     x = Sx[i] = xx + HCompl->Julia_Parm[X];
  591.     y = Sy[i] = yy + HCompl->Julia_Parm[Y];
  592.     z = Sz[i] = zz + HCompl->Julia_Parm[Z];
  593.     w = Sw[i] = ww + HCompl->Julia_Parm[T];
  594.   }
  595.  
  596.   *Dist = Precision;
  597.  
  598.   return (TRUE);
  599. }
  600.  
  601.  
  602.  
  603. /*****************************************************************************
  604. *
  605. * FUNCTION
  606. *
  607. * INPUT
  608. *
  609. * OUTPUT
  610. *
  611. * RETURNS
  612. *
  613. * AUTHOR
  614. *
  615. *   Pascal Massimino
  616. *
  617. * DESCRIPTION
  618. *
  619. * CHANGES
  620. *
  621. ******************************************************************************/
  622.  
  623. void Normal_Calc_HCompl_z3(Result, N_Max, fractal)
  624. VECTOR Result;
  625. int N_Max;
  626. FRACTAL *fractal;
  627. {
  628.   DBL n1, n2, n3, n4;
  629.   int i;
  630.   DBL x, y, z, w;
  631.   DBL xx, yy, zz, ww;
  632.  
  633.   /*
  634.    * Algebraic properties of hypercomplexes allows simplifications in
  635.    * computations...
  636.    */
  637.  
  638.   x = Sx[0];
  639.   y = Sy[0];
  640.   z = Sz[0];
  641.   w = Sw[0];
  642.  
  643.   for (i = 1; i < N_Max; ++i)
  644.   {
  645.     /*
  646.      * For a map z->f(z), f depending on c, one must perform here :
  647.      * 
  648.      * (x,y,z,w) * df/dz(Sx[i],Sy[i],Sz[i],Sw[i]) -> (x,y,z,w) ,
  649.      * 
  650.      * up to a constant.
  651.      */
  652.  
  653.     /******************* Case z->z^2+c *****************/
  654.  
  655.     /* the df/dz part needs no work */
  656.  
  657.     HMult(xx, yy, zz, ww, Sx[i], Sy[i], Sz[i], Sw[i], x, y, z, w);
  658.  
  659.     x = xx;
  660.     y = yy;
  661.     z = zz;
  662.     w = ww;
  663.   }
  664.  
  665.   n1 = Sx[N_Max];
  666.   n2 = Sy[N_Max];
  667.   n3 = Sz[N_Max];
  668.   n4 = Sw[N_Max];
  669.  
  670.   Result[X] = x * n1 + y * n2 + z * n3 + w * n4;
  671.   Result[Y] = -y * n1 + x * n2 - w * n3 + z * n4;
  672.   Result[Z] = -z * n1 - w * n2 + x * n3 + y * n4;
  673. }
  674.  
  675.  
  676.  
  677. /*****************************************************************************
  678. *
  679. * FUNCTION
  680. *
  681. * INPUT
  682. *
  683. * OUTPUT
  684. *
  685. * RETURNS
  686. *
  687. * AUTHOR
  688. *
  689. *   Pascal Massimino
  690. *
  691. * DESCRIPTION
  692. *
  693. * CHANGES
  694. *
  695. ******************************************************************************/
  696.  
  697. int F_Bound_HCompl_z3(Ray, Fractal, Depth_Min, Depth_Max)
  698. RAY *Ray;
  699. FRACTAL *Fractal;
  700. DBL *Depth_Min;
  701. DBL *Depth_Max;
  702. {
  703.   return F_Bound_HCompl(Ray, Fractal, Depth_Min, Depth_Max);
  704. }
  705.  
  706. /*--------------------------- Inv -------------------------------*/
  707.  
  708.  
  709. /*****************************************************************************
  710. *
  711. * FUNCTION
  712. *
  713. * INPUT
  714. *
  715. * OUTPUT
  716. *
  717. * RETURNS
  718. *
  719. * AUTHOR
  720. *
  721. *   Pascal Massimino
  722. *
  723. * DESCRIPTION
  724. *
  725. * CHANGES
  726. *
  727. ******************************************************************************/
  728.  
  729. int Iteration_HCompl_Reciprocal(IPoint, HCompl)
  730. VECTOR IPoint;
  731. FRACTAL *HCompl;
  732. {
  733.   int i;
  734.   DBL Norm, xx, yy, zz, ww;
  735.   DBL Exit_Value;
  736.   DBL x, y, z, w;
  737.  
  738.   x = Sx[0] = IPoint[X];
  739.   y = Sy[0] = IPoint[Y];
  740.   z = Sz[0] = IPoint[Z];
  741.   w = Sw[0] = (HCompl->SliceDist
  742.                   - HCompl->Slice[X]*x 
  743.                   - HCompl->Slice[Y]*y 
  744.                   - HCompl->Slice[Z]*z)/HCompl->Slice[T]; 
  745.  
  746.   Exit_Value = HCompl->Exit_Value;
  747.  
  748.   for (i = 1; i <= HCompl->n; ++i)
  749.   {
  750.     Norm = x * x + y * y + z * z + w * w;
  751.  
  752.     if (Norm > Exit_Value)
  753.     {
  754.       return (FALSE);
  755.     }
  756.  
  757.     HReciprocal(&xx, &yy, &zz, &ww, x, y, z, w);
  758.  
  759.     x = Sx[i] = xx + HCompl->Julia_Parm[X];
  760.     y = Sy[i] = yy + HCompl->Julia_Parm[Y];
  761.     z = Sz[i] = zz + HCompl->Julia_Parm[Z];
  762.     w = Sw[i] = ww + HCompl->Julia_Parm[T];
  763.  
  764.   }
  765.  
  766.   return (TRUE);
  767. }
  768.  
  769.  
  770.  
  771. /*****************************************************************************
  772. *
  773. * FUNCTION
  774. *
  775. * INPUT
  776. *
  777. * OUTPUT
  778. *
  779. * RETURNS
  780. *
  781. * AUTHOR
  782. *
  783. *   Pascal Massimino
  784. *
  785. * DESCRIPTION
  786. *
  787. * CHANGES
  788. *
  789. ******************************************************************************/
  790.  
  791. int D_Iteration_HCompl_Reciprocal(IPoint, HCompl, Dist)
  792. VECTOR IPoint;
  793. FRACTAL *HCompl;
  794. DBL *Dist;
  795. {
  796.   int i;
  797.   DBL xx, yy, zz, ww;
  798.   DBL Exit_Value, F_Value, Step;
  799.   DBL x, y, z, w;
  800.   VECTOR H_Normal;
  801.  
  802.   x = Sx[0] = IPoint[X];
  803.   y = Sy[0] = IPoint[Y];
  804.   z = Sz[0] = IPoint[Z];
  805.   w = Sw[0] = (HCompl->SliceDist
  806.                   - HCompl->Slice[X]*x 
  807.                   - HCompl->Slice[Y]*y 
  808.                   - HCompl->Slice[Z]*z)/HCompl->Slice[T]; 
  809.  
  810.   Exit_Value = HCompl->Exit_Value;
  811.  
  812.   for (i = 1; i <= HCompl->n; ++i)
  813.   {
  814.     F_Value = x * x + y * y + z * z + w * w;
  815.  
  816.     if (F_Value > Exit_Value)
  817.     {
  818.       Normal_Calc_HCompl_Reciprocal(H_Normal, i - 1, HCompl);
  819.  
  820.       VDot(Step, H_Normal, Direction);
  821.  
  822.       if (Step < -Fractal_Tolerance)
  823.       {
  824.         Step = -2.0 * Step;
  825.  
  826.         if ((F_Value > Precision * Step) && F_Value < (30 * Precision * Step))
  827.         {
  828.           *Dist = F_Value / Step;
  829.  
  830.           return (FALSE);
  831.         }
  832.       }
  833.  
  834.       *Dist = Precision;
  835.  
  836.       return (FALSE);
  837.     }
  838.  
  839.     HReciprocal(&xx, &yy, &zz, &ww, x, y, z, w);
  840.  
  841.     x = Sx[i] = xx + HCompl->Julia_Parm[X];
  842.     y = Sy[i] = yy + HCompl->Julia_Parm[Y];
  843.     z = Sz[i] = zz + HCompl->Julia_Parm[Z];
  844.     w = Sw[i] = ww + HCompl->Julia_Parm[T];
  845.  
  846.   }
  847.  
  848.   *Dist = Precision;
  849.  
  850.   return (TRUE);
  851. }
  852.  
  853.  
  854.  
  855. /*****************************************************************************
  856. *
  857. * FUNCTION
  858. *
  859. * INPUT
  860. *
  861. * OUTPUT
  862. *
  863. * RETURNS
  864. *
  865. * AUTHOR
  866. *
  867. *   Pascal Massimino
  868. *
  869. * DESCRIPTION
  870. *
  871. * CHANGES
  872. *
  873. ******************************************************************************/
  874.  
  875. void Normal_Calc_HCompl_Reciprocal(Result, N_Max, fractal)
  876. VECTOR Result;
  877. int N_Max;
  878. FRACTAL *fractal;
  879. {
  880.   DBL n1, n2, n3, n4;
  881.   int i;
  882.   DBL x, y, z, w;
  883.   DBL xx, yy, zz, ww;
  884.   DBL xxx, yyy, zzz, www;
  885.  
  886.   /*
  887.    * Algebraic properties of hypercomplexes allows simplifications in
  888.    * computations...
  889.    */
  890.  
  891.   x = Sx[0];
  892.   y = Sy[0];
  893.   z = Sz[0];
  894.   w = Sw[0];
  895.  
  896.   for (i = 1; i < N_Max; ++i)
  897.   {
  898.     /******************* Case: z->1/z+c *****************/
  899.  
  900.     HReciprocal(&xx, &yy, &zz, &ww, Sx[i], Sy[i], Sz[i], Sw[i]);
  901.  
  902.     HSqr(xxx, yyy, zzz, www, xx, yy, zz, ww);
  903.  
  904.     HMult(xx, yy, zz, ww, x, y, z, w, -xxx, -yyy, -zzz, -www);
  905.  
  906.     x = xx;
  907.     y = yy;
  908.     z = zz;
  909.     w = ww;
  910.   }
  911.  
  912.   n1 = Sx[N_Max];
  913.   n2 = Sy[N_Max];
  914.   n3 = Sz[N_Max];
  915.   n4 = Sw[N_Max];
  916.  
  917.   Result[X] = x * n1 + y * n2 + z * n3 + w * n4;
  918.   Result[Y] = -y * n1 + x * n2 - w * n3 + z * n4;
  919.   Result[Z] = -z * n1 - w * n2 + x * n3 + y * n4;
  920. }
  921.  
  922.  
  923.  
  924. /*****************************************************************************
  925. *
  926. * FUNCTION
  927. *
  928. * INPUT
  929. *
  930. * OUTPUT
  931. *
  932. * RETURNS
  933. *
  934. * AUTHOR
  935. *
  936. *   Pascal Massimino
  937. *
  938. * DESCRIPTION
  939. *
  940. * CHANGES
  941. *
  942. ******************************************************************************/
  943.  
  944. int F_Bound_HCompl_Reciprocal(Ray, Fractal, Depth_Min, Depth_Max)
  945. RAY *Ray;
  946. FRACTAL *Fractal;
  947. DBL *Depth_Min, *Depth_Max;
  948. {
  949.   return F_Bound_HCompl(Ray, Fractal, Depth_Min, Depth_Max);
  950. }
  951.  
  952. /*--------------------------- Func -------------------------------*/
  953.  
  954.  
  955. /*****************************************************************************
  956. *
  957. * FUNCTION
  958. *
  959. * INPUT
  960. *
  961. * OUTPUT
  962. *
  963. * RETURNS
  964. *
  965. * AUTHOR
  966. *
  967. *   Pascal Massimino
  968. *
  969. * DESCRIPTION
  970. *
  971. * CHANGES
  972. *
  973. ******************************************************************************/
  974.  
  975. int Iteration_HCompl_Func(IPoint, HCompl)
  976. VECTOR IPoint;
  977. FRACTAL *HCompl;
  978. {
  979.   int i;
  980.   DBL Norm, xx, yy, zz, ww;
  981.   DBL Exit_Value;
  982.   DBL x, y, z, w;
  983.  
  984.   x = Sx[0] = IPoint[X];
  985.   y = Sy[0] = IPoint[Y];
  986.   z = Sz[0] = IPoint[Y];
  987.   w = Sw[0] = (HCompl->SliceDist
  988.                   - HCompl->Slice[X]*x 
  989.                   - HCompl->Slice[Y]*y 
  990.                   - HCompl->Slice[Z]*z)/HCompl->Slice[T]; 
  991.  
  992.   Exit_Value = HCompl->Exit_Value;
  993.  
  994.   for (i = 1; i <= HCompl->n; ++i)
  995.   {
  996.     Norm = x * x + y * y + z * z + w * w;
  997.  
  998.     if (Norm > Exit_Value)
  999.     {
  1000.       return (FALSE);
  1001.     }
  1002.  
  1003.     HFunc(&xx, &yy, &zz, &ww, x, y, z, w, HCompl);
  1004.  
  1005.     x = Sx[i] = xx + HCompl->Julia_Parm[X];
  1006.     y = Sy[i] = yy + HCompl->Julia_Parm[Y];
  1007.     z = Sz[i] = zz + HCompl->Julia_Parm[Z];
  1008.     w = Sw[i] = ww + HCompl->Julia_Parm[T];
  1009.  
  1010.   }
  1011.  
  1012.   return (TRUE);
  1013. }
  1014.  
  1015.  
  1016.  
  1017. /*****************************************************************************
  1018. *
  1019. * FUNCTION
  1020. *
  1021. * INPUT
  1022. *
  1023. * OUTPUT
  1024. *
  1025. * RETURNS
  1026. *
  1027. * AUTHOR
  1028. *
  1029. *   Pascal Massimino
  1030. *
  1031. * DESCRIPTION
  1032. *
  1033. * CHANGES
  1034. *
  1035. ******************************************************************************/
  1036.  
  1037. int D_Iteration_HCompl_Func(IPoint, HCompl, Dist)
  1038. VECTOR IPoint;
  1039. FRACTAL *HCompl;
  1040. DBL *Dist;
  1041. {
  1042.   int i;
  1043.   DBL xx, yy, zz, ww;
  1044.   DBL Exit_Value, F_Value, Step;
  1045.   DBL x, y, z, w;
  1046.   VECTOR H_Normal;
  1047.  
  1048.   x = Sx[0] = IPoint[X];
  1049.   y = Sy[0] = IPoint[Y];
  1050.   z = Sz[0] = IPoint[Z];
  1051.   w = Sw[0] = (HCompl->SliceDist
  1052.                   - HCompl->Slice[X]*x 
  1053.                   - HCompl->Slice[Y]*y 
  1054.                   - HCompl->Slice[Z]*z)/HCompl->Slice[T]; 
  1055.  
  1056.   Exit_Value = HCompl->Exit_Value;
  1057.  
  1058.   for (i = 1; i <= HCompl->n; ++i)
  1059.   {
  1060.     F_Value = x * x + y * y + z * z + w * w;
  1061.  
  1062.     if (F_Value > Exit_Value)
  1063.     {
  1064.       Normal_Calc_HCompl_Func(H_Normal, i - 1, HCompl);
  1065.  
  1066.       VDot(Step, H_Normal, Direction);
  1067.  
  1068.       if (Step < -Fractal_Tolerance)
  1069.       {
  1070.         Step = -2.0 * Step;
  1071.  
  1072.         if ((F_Value > Precision * Step) && F_Value < (30 * Precision * Step))
  1073.         {
  1074.           *Dist = F_Value / Step;
  1075.  
  1076.           return (FALSE);
  1077.         }
  1078.       }
  1079.  
  1080.       *Dist = Precision;
  1081.  
  1082.       return (FALSE);
  1083.     }
  1084.  
  1085.     HFunc(&xx, &yy, &zz, &ww, x, y, z, w, HCompl);
  1086.  
  1087.     x = Sx[i] = xx + HCompl->Julia_Parm[X];
  1088.     y = Sy[i] = yy + HCompl->Julia_Parm[Y];
  1089.     z = Sz[i] = zz + HCompl->Julia_Parm[Z];
  1090.     w = Sw[i] = ww + HCompl->Julia_Parm[T];
  1091.  
  1092.   }
  1093.  
  1094.   *Dist = Precision;
  1095.  
  1096.   return (TRUE);
  1097. }
  1098.  
  1099.  
  1100.  
  1101. /*****************************************************************************
  1102. *
  1103. * FUNCTION
  1104. *
  1105. * INPUT
  1106. *
  1107. * OUTPUT
  1108. *
  1109. * RETURNS
  1110. *
  1111. * AUTHOR
  1112. *
  1113. *   Pascal Massimino
  1114. *
  1115. * DESCRIPTION
  1116. *
  1117. * CHANGES
  1118. *
  1119. ******************************************************************************/
  1120.  
  1121. void Normal_Calc_HCompl_Func(Result, N_Max, fractal)
  1122. VECTOR Result;
  1123. int N_Max;
  1124. FRACTAL *fractal;
  1125. {
  1126.   DBL n1, n2, n3, n4;
  1127.   int i;
  1128.   DBL x, y, z, w;
  1129.   DBL xx, yy, zz, ww;
  1130.   DBL xxx, yyy, zzz, www;
  1131.  
  1132.   /*
  1133.    * Algebraic properties of hypercomplexes allows simplifications in
  1134.    * computations...
  1135.    */
  1136.  
  1137.   x = Sx[0];
  1138.   y = Sy[0];
  1139.   z = Sz[0];
  1140.   w = Sw[0];
  1141.  
  1142.   for (i = 1; i < N_Max; ++i)
  1143.   {
  1144.     /**************** Case: z-> f(z)+c ************************/
  1145.  
  1146.     HFunc(&xx, &yy, &zz, &ww, Sx[i], Sy[i], Sz[i], Sw[i], fractal);
  1147.  
  1148.     HMult(xxx, yyy, zzz, www, xx, yy, zz, ww, x, y, z, w);
  1149.  
  1150.     x = xxx;
  1151.     y = yyy;
  1152.     z = zzz;
  1153.     w = www;
  1154.   }
  1155.  
  1156.   n1 = Sx[N_Max];
  1157.   n2 = Sy[N_Max];
  1158.   n3 = Sz[N_Max];
  1159.   n4 = Sw[N_Max];
  1160.  
  1161.   Result[X] = x * n1 + y * n2 + z * n3 + w * n4;
  1162.   Result[Y] = -y * n1 + x * n2 - w * n3 + z * n4;
  1163.   Result[Z] = -z * n1 - w * n2 + x * n3 + y * n4;
  1164. }
  1165.  
  1166.  
  1167.  
  1168. /*****************************************************************************
  1169. *
  1170. * FUNCTION
  1171. *
  1172. * INPUT
  1173. *
  1174. * OUTPUT
  1175. *
  1176. * RETURNS
  1177. *
  1178. * AUTHOR
  1179. *
  1180. *   Pascal Massimino
  1181. *
  1182. * DESCRIPTION
  1183. *
  1184. * CHANGES
  1185. *
  1186. ******************************************************************************/
  1187.  
  1188. int F_Bound_HCompl_Func(Ray, Fractal, Depth_Min, Depth_Max)
  1189. RAY *Ray;
  1190. FRACTAL *Fractal;
  1191. DBL *Depth_Min, *Depth_Max;
  1192. {
  1193.   return F_Bound_HCompl(Ray, Fractal, Depth_Min, Depth_Max);
  1194. }
  1195.  
  1196. /*****************************************************************************
  1197. *
  1198. * FUNCTION  Complex transcental functions
  1199. *
  1200. * INPUT     pointer to source complex number
  1201. *
  1202. * OUTPUT    fn(input)
  1203. *
  1204. * RETURNS   void
  1205. *
  1206. * AUTHOR
  1207. *
  1208. *   Tim Wegner
  1209. *
  1210. * DESCRIPTION  Calculate common functions on complexes
  1211. *   Since our purpose is fractals, error checking is lax.            
  1212. *
  1213. * CHANGES
  1214. *
  1215. ******************************************************************************/
  1216.  
  1217. void Complex_Mult (target, source1, source2)
  1218. CMPLX *target;
  1219. CMPLX *source1;
  1220. CMPLX *source2;
  1221. {
  1222.   DBL tmpx;
  1223.   tmpx = source1->x * source2->x - source1->y * source2->y;
  1224.   target->y = source1->x * source2->y + source1->y * source2->x;
  1225.   target->x = tmpx;
  1226. }
  1227.  
  1228. void Complex_Div (target, source1, source2)
  1229. CMPLX *target;
  1230. CMPLX *source1;
  1231. CMPLX *source2;
  1232. {
  1233.   DBL mod,tmpx,yxmod,yymod;
  1234.   mod = Sqr(source2->x) + Sqr(source2->y);
  1235.   if (mod==0)
  1236.      return;
  1237.   yxmod = source2->x/mod;
  1238.   yymod = - source2->y/mod;
  1239.   tmpx = source1->x * yxmod - source1->y * yymod;
  1240.   target->y = source1->x * yymod + source1->y * yxmod;
  1241.   target->x = tmpx;
  1242. } /* End Complex_Mult() */
  1243.  
  1244. void Complex_Exp (target, source)
  1245. CMPLX *target;
  1246. CMPLX *source;
  1247. {
  1248.   DBL expx;
  1249.   expx = exp(source->x);
  1250.   target->x = expx * cos(source->y);
  1251.   target->y = expx * sin(source->y);
  1252. } /* End Complex_Exp() */
  1253.  
  1254. void Complex_Sin (target, source)
  1255. CMPLX *target;
  1256. CMPLX *source;
  1257. {
  1258.   target->x = sin(source->x) * cosh(source->y);
  1259.   target->y = cos(source->x) * sinh(source->y);
  1260. } /* End Complex_Sin() */
  1261.  
  1262. void Complex_Sinh (target, source)
  1263. CMPLX *target;
  1264. CMPLX *source;
  1265. {
  1266.   target->x = sinh(source->x) * cos(source->y);
  1267.   target->y = cosh(source->x) * sin(source->y);
  1268. } /* End Complex_Sinh() */
  1269.  
  1270.  
  1271. void Complex_Cos (target, source)
  1272. CMPLX *target;
  1273. CMPLX *source;
  1274. {
  1275.   target->x = cos(source->x) * cosh(source->y);
  1276.   target->y = -sin(source->x) * sinh(source->y);
  1277. } /* End Complex_Cos() */
  1278.  
  1279. void Complex_Cosh (target, source)
  1280. CMPLX *target;
  1281. CMPLX *source;
  1282. {
  1283.   target->x = cosh(source->x) * cos(source->y);
  1284.   target->y = sinh(source->x) * sin(source->y);
  1285. } /* End Complex_Cosh() */
  1286.  
  1287.  
  1288. void Complex_Log (target, source)
  1289. CMPLX *target;
  1290. CMPLX *source;
  1291. {
  1292.   DBL mod,zx,zy;
  1293.   mod = sqrt(source->x * source->x + source->y * source->y);
  1294.   zx = log(mod);
  1295.   zy = atan2(source->y,source->x);
  1296.  
  1297.   target->x = zx;
  1298.   target->y = zy;
  1299. } /* End Complex_Log() */
  1300.  
  1301. void Complex_Sqrt(target, source)
  1302. CMPLX *target;
  1303. CMPLX *source;
  1304. {
  1305.   DBL mag;
  1306.   DBL theta;
  1307.  
  1308.   if(source->x == 0.0 && source->y == 0.0)
  1309.   {
  1310.     target->x = target->y = 0.0;
  1311.   }
  1312.   else
  1313.   {   
  1314.     mag   = sqrt(sqrt(Sqr(source->x) + Sqr(source->y)));
  1315.     theta = atan2(source->y, source->x) / 2;
  1316.     target->y = mag * sin(theta);
  1317.     target->x = mag * cos(theta);
  1318.   }
  1319. } /* End Complex_Sqrt() */
  1320.  
  1321. /* rz=Arcsin(z)=-i*Log{i*z+sqrt(1-z*z)} */
  1322. void Complex_ASin(target, source)
  1323. CMPLX *target;
  1324. CMPLX *source;
  1325. {
  1326.   CMPLX tempz1,tempz2;
  1327.  
  1328.   Complex_Mult(&tempz1, source, source);
  1329.   tempz1.x = 1 - tempz1.x; tempz1.y = -tempz1.y; 
  1330.   Complex_Sqrt( &tempz1, &tempz1);
  1331.  
  1332.   tempz2.x = -source->y; tempz2.y = source->x;
  1333.   tempz1.x += tempz2.x;  tempz1.y += tempz2.y;
  1334.  
  1335.   Complex_Log( &tempz1, &tempz1);
  1336.   target->x = tempz1.y;  target->y = -tempz1.x;   
  1337. }   /* End Complex_ASin() */
  1338.  
  1339.  
  1340. void Complex_ACos(target, source)
  1341. CMPLX *target;
  1342. CMPLX *source;
  1343. {
  1344.   CMPLX temp;
  1345.  
  1346.   Complex_Mult(&temp, source, source);
  1347.   temp.x -= 1;
  1348.   Complex_Sqrt(&temp, &temp);
  1349.  
  1350.   temp.x += source->x; temp.y += source->y;
  1351.   
  1352.   Complex_Log(&temp, &temp);
  1353.   target->x = temp.y;  target->y = -temp.x; 
  1354. }   /* End Complex_ACos() */
  1355.  
  1356. void Complex_ASinh(target, source)
  1357. CMPLX *target;
  1358. CMPLX *source;
  1359. {
  1360.   CMPLX temp;
  1361.  
  1362.   Complex_Mult (&temp, source, source);
  1363.   temp.x += 1;   
  1364.   Complex_Sqrt (&temp, &temp);
  1365.   temp.x += source->x; temp.y += source->y; 
  1366.   Complex_Log(target, &temp);
  1367. }  /* End Complex_ASinh */
  1368.  
  1369. /* rz=Arccosh(z)=Log(z+sqrt(z*z-1)} */
  1370. void Complex_ACosh (target, source)
  1371. CMPLX *target;
  1372. CMPLX *source;
  1373. {
  1374.   CMPLX tempz;
  1375.   Complex_Mult(&tempz, source, source);
  1376.   tempz.x -= 1; 
  1377.   Complex_Sqrt (&tempz, &tempz);
  1378.   tempz.x = source->x + tempz.x; tempz.y = source->y + tempz.y;
  1379.   Complex_Log (target, &tempz);
  1380. }   /* End Complex_ACosh() */
  1381.  
  1382. /* rz=Arctanh(z)=1/2*Log{(1+z)/(1-z)} */
  1383. void Complex_ATanh(target, source)
  1384. CMPLX *target;
  1385. CMPLX *source;
  1386. {
  1387.   CMPLX temp0,temp1,temp2;
  1388.    
  1389.   if( source->x == 0.0)
  1390.   {
  1391.     target->x = 0;
  1392.     target->y = atan( source->y);
  1393.     return;
  1394.   }
  1395.   else
  1396.   {
  1397.     if( fabs(source->x) == 1.0 && source->y == 0.0)
  1398.     {
  1399.       return;
  1400.     }
  1401.     else if( fabs( source->x) < 1.0 && source->y == 0.0)
  1402.     {
  1403.       target->x = log((1+source->x)/(1-source->x))/2;
  1404.       target->y = 0;
  1405.       return;
  1406.     } 
  1407.     else
  1408.     {
  1409.       temp0.x = 1 + source->x; temp0.y = source->y;
  1410.       temp1.x = 1 - source->x; temp1.y = -source->y; 
  1411.       Complex_Div(&temp2, &temp0, &temp1);
  1412.       Complex_Log(&temp2, &temp2);
  1413.       target->x = .5 * temp2.x; target->y = .5 * temp2.y;
  1414.       return;
  1415.     }
  1416.   }
  1417. }   /* End Complex_ATanh() */
  1418.  
  1419. /* rz=Arctan(z)=i/2*Log{(1-i*z)/(1+i*z)} */
  1420. void Complex_ATan(target, source)
  1421. CMPLX *target;
  1422. CMPLX *source;
  1423. {
  1424.   CMPLX temp0,temp1,temp2,temp3;
  1425.   if( source->x == 0.0 && source->y == 0.0)
  1426.     target->x = target->y = 0;
  1427.   else if( source->x != 0.0 && source->y == 0.0){
  1428.     target->x = atan(source->x);
  1429.     target->y = 0;
  1430.   }
  1431.   else if( source->x == 0.0 && source->y != 0.0){
  1432.     temp0.x = source->y;  temp0.y = 0.0;
  1433.     Complex_ATanh(&temp0, &temp0);
  1434.     target->x = -temp0.y; target->y = temp0.x;
  1435.   } 
  1436.   else if( source->x != 0.0 && source->y != 0.0)
  1437.   {
  1438.     temp0.x = -source->y; temp0.y = source->x; 
  1439.     temp1.x = 1 - temp0.x; temp1.y = -temp0.y;   
  1440.     temp2.x = 1 + temp0.x; temp2.y = temp0.y; 
  1441.  
  1442.     Complex_Div(&temp3, &temp1, &temp2);
  1443.     Complex_Log(&temp3, &temp3);
  1444.     target->x = -temp3.y * .5; target->y = .5 * temp3.x; 
  1445.   }
  1446. }   /* End Complex_ATanz() */
  1447.  
  1448. void Complex_Tan(target, source)
  1449. CMPLX *target;
  1450. CMPLX *source;
  1451. {
  1452.   DBL x, y, sinx, cosx, sinhy, coshy, denom;
  1453.   x = 2 * source->x;
  1454.   y = 2 * source->y;
  1455.   sinx = sin(x); cosx = cos(x);
  1456.   sinhy = sinh(y); coshy = cosh(y);
  1457.   denom = cosx + coshy;
  1458.   if(denom == 0)
  1459.     return;
  1460.   target->x = sinx/denom;
  1461.   target->y = sinhy/denom;
  1462. }   /* End Complex_Tan() */
  1463.  
  1464. void Complex_Tanh(target, source)
  1465. CMPLX *target;
  1466. CMPLX *source;
  1467. {
  1468.   DBL x, y, siny, cosy, sinhx, coshx, denom;
  1469.   x = 2 * source->x;
  1470.   y = 2 * source->y;
  1471.   siny = sin(y); cosy = cos(y);
  1472.   sinhx = sinh(x); coshx = cosh(x);
  1473.   denom = coshx + cosy;
  1474.   if(denom == 0)
  1475.     return;
  1476.   target->x = sinhx/denom;
  1477.   target->y = siny/denom;
  1478. }   /* End Complex_Tanh() */
  1479.  
  1480.  
  1481. void Complex_Power(target, source1, source2) 
  1482. CMPLX *target;
  1483. CMPLX *source1;
  1484. CMPLX *source2;
  1485. {
  1486.   CMPLX cLog, t;
  1487.   DBL e2x;
  1488.  
  1489.   if(source1->x == 0 && source1->y == 0) 
  1490.   {
  1491.     target->x = target->y = 0.0;
  1492.     return;
  1493.   }
  1494.  
  1495.   Complex_Log (&cLog, source1);
  1496.   Complex_Mult (&t, &cLog, source2);
  1497.  
  1498.   if(t.x < -690)
  1499.     e2x = 0;
  1500.   else
  1501.     e2x = exp(t.x);
  1502.   target->x = e2x * cos(t.y);
  1503.   target->y = e2x * sin(t.y);
  1504. }
  1505.  
  1506. #if 1
  1507. void Complex_Pwr(target, source) 
  1508. CMPLX *target;
  1509. CMPLX *source;
  1510. {
  1511.   Complex_Power(target, source, &exponent);
  1512. }
  1513.  
  1514. #else
  1515. void Complex_Pwr(target, source) 
  1516. CMPLX *target;
  1517. CMPLX *source;
  1518. {
  1519.   /* the sqr dunction for testing */
  1520.   Complex_Mult(target, source, source);
  1521. }
  1522. #endif
  1523.